openHDF5
Megan Cattau
June 20, 2016
# load the required packages
library(raster)## Loading required package: splibrary(rhdf5)
library(rgdal)## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3
We’re going to open the reflectance dataset from the Teakettle site and look at the structure
# if you type "../" then hit tab, the file drop-down options will open automatically
f<-"../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"
# view H5 structure to explore it
h5ls(f)## group name otype dclass
## 0 / ATCOR_Input_File H5I_DATASET STRING
## 1 / ATCOR_Processing_Log H5I_DATASET STRING
## 2 / Aerosol Optical Depth H5I_DATASET INTEGER
## 3 / Aspect H5I_DATASET FLOAT
## 4 / Cast Shadow H5I_DATASET INTEGER
## 5 / Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6 / Haze-Cloud-Water Map H5I_DATASET INTEGER
## 7 / Illumination Factor H5I_DATASET INTEGER
## 8 / Path Length H5I_DATASET FLOAT
## 9 / Processing Version H5I_DATASET STRING
## 10 / Reflectance H5I_DATASET INTEGER
## 11 / Shadow_Processing_Log H5I_DATASET STRING
## 12 / Sky View Factor H5I_DATASET INTEGER
## 13 / Skyview_Processing_Log H5I_DATASET STRING
## 14 / Slope H5I_DATASET FLOAT
## 15 / Slope_Aspect_Processing_Log H5I_DATASET STRING
## 16 / Solar Azimuth Angle H5I_DATASET FLOAT
## 17 / Solar Zenith Angle H5I_DATASET FLOAT
## 18 / Surface Elevation H5I_DATASET FLOAT
## 19 / Visibility Index Map H5I_DATASET INTEGER
## 20 / Water Vapor Column H5I_DATASET INTEGER
## 21 / coordinate system string H5I_DATASET STRING
## 22 / flightAltitude H5I_DATASET FLOAT
## 23 / flightHeading H5I_DATASET FLOAT
## 24 / flightTime H5I_DATASET FLOAT
## 25 / fwhm H5I_DATASET FLOAT
## 26 / map info H5I_DATASET STRING
## 27 / to-sensor azimuth angle H5I_DATASET FLOAT
## 28 / to-sensor zenith angle H5I_DATASET FLOAT
## 29 / wavelength H5I_DATASET FLOAT
## dim
## 0 1
## 1 1
## 2 544 x 578 x 1
## 3 544 x 578 x 1
## 4 544 x 578 x 1
## 5 544 x 578 x 1
## 6 544 x 578 x 1
## 7 544 x 578 x 1
## 8 544 x 578 x 1
## 9 1
## 10 544 x 578 x 426
## 11 1
## 12 544 x 578 x 1
## 13 1
## 14 544 x 578 x 1
## 15 1
## 16 1 x 1
## 17 1 x 1
## 18 544 x 578 x 1
## 19 544 x 578 x 1
## 20 544 x 578 x 1
## 21 1
## 22 5732053
## 23 5732053
## 24 5732053
## 25 426 x 1
## 26 1
## 27 544 x 578 x 1
## 28 544 x 578 x 1
## 29 426 x 1
Import spatial information - view the attributes (metadata) for the map info dataset within the NEON H5 file.
# tell R where the data are located
# in HDF viewer, the "map info" group tells you the coordinate reference system and the coordinate loaction of the upper left corner of the file and spatial resolution
mapInfo<-h5read(f,
"map info", # one can find the name of this data in HDF view
read.attributes = TRUE)
mapInfo## [1] "UTM,1,1,325963.0,4103482.0,1.0000000000e+000,1.0000000000e+000,11,North,WGS-84,units=Meters"
## attr(,"Description")
## [1] "Basic Map information for envi style programs"
We’re going to import the reflectance metadata and define the scale factor and no data value
# read in reflectance data attributes
reflInfo<-h5readAttributes(f,
name = "Reflectance")
reflInfo## $DIMENSION_LABELS
## [1] "Wavelength" "Line" "Sample"
##
## $Description
## [1] "Atmospherically corrected reflectance."
##
## $`Scale Factor`
## [1] 10000
##
## $Unit
## [1] "unitless. Valid range 0-1."
##
## $`data ignore value`
## [1] "15000.0"# define scale factor and no data value as objects
scaleFactor<-reflInfo$`Scale Factor`
noDataValue<-as.numeric(reflInfo$`data ignore value`)
Import data dimensions
# open file for viewing. This is a direct connection to the file, but it doesn't read it in. We don't want to read in a huge data cube into memory.
fid<-H5Fopen(f)
# Use this to open the reflectance dataset itself
did<-H5Dopen(fid, "Reflectance")
did## HDF5 DATASET
## name /Reflectance
## filename
## type H5T_STD_I16LE
## rank 3
## size 544 x 578 x 426
## maxsize 544 x 578 x 426# this has columns first rather than rows
# grab the dimensions of the object
# you could also extract dimensions from H5ls, but this creates a text object and then you have to split it out
# H5ls useful if you want to start looping through multiple files
sid <- H5Dget_space(did)
dims <- H5Sget_simple_extent_dims(sid)$size
dims## [1] 544 578 426# Columns are the FIRST dimension, then rows. wavelength is the THIRD dimension rather than the first.
# close all open connections
H5Sclose(sid)
H5Dclose(did)
H5Fclose(fid)
Read Reflectance Data
Once we know the dimensions of the data, we can more efficiently slice out chunks or subsets of our data out. The power of HDF5 is that it allows us to store large heterogeneous data. However, we can quickly and efficiently access those data through “slicing” or extracting subsets of the data.
# Extract or "slice" data (one layer of datacube) for band 56 from the HDF5 file
b56<- h5read(f,"Reflectance", index=list(1:dims[1],1:dims[2],56))
# note the data come in as an array
class(b56)## [1] "array"
Convert to Matrix and plot the data
# Convert from array to matrix so we can plot and convert to a raster
b56 <- b56[,,1]
# plot the data
image(b56)
# stretch it
image(log(b56),
main="Band 56 with log Transformation")
# view distribution of reflectance data
# force non scientific notation
options("scipen"=100, "digits"=4)
hist(b56,
col="springgreen",
main="Distribution of Reflectance Values \nBand 56")
Data Clean-up Set the no data value to 15,000. Apply the scale factor to the data (10,000).
# extract no data value from the attributes
noDataVal <- as.integer(reflInfo$`data ignore value`)
# set all reflectance values = 15,000 to NA
b56[b56 == noDataVal] <- NA
# Extract the scale factor as an object
scaleFactor <- reflInfo$`Scale Factor`
# divide all values in our B56 object by the scale factor to get a range of
# reflectance values between 0-1 (the valid range)
b56 <- b56/scaleFactor
# view distribution of reflectance values
hist(b56,
main="Distribution with NoData Value Considered\nData Scaled")
Unflip the data
b56<-t(b56)
image(log(b56), main="Band 56\nTransposed Values")
Create Spatial Extent
mapInfo<-unlist(strsplit(mapInfo, ","))
# X,Y left corner coordinate - change to numeric
xMin <- as.numeric(mapInfo[4])
yMax <- as.numeric(mapInfo[5])
# get x and y resolution
xres <- as.numeric(mapInfo[6])
yres <- as.numeric(mapInfo[7])
# finally calculate the xMax value and the yMin value from the dimensions
# we grabbed above. The xMax is the left corner + number of columns* resolution.
xMax <- xMin + (dims[1]*xres)
yMin <- yMax - (dims[2]*yres)
# Define the raster extent
rasExt <- extent(xMin, xMax,yMin,yMax)
# Create a raster and assign it it's spatial extent
b56r <- raster(b56,
crs=CRS("+init=epsg:32611"))
# assign CRS
extent(b56r) <- rasExt
# view raster object attributes
b56r## class : RasterLayer
## dimensions : 578, 544, 314432 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 325963, 326507, 4102904, 4103482 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 0.0033, 0.5678 (min, max)# plot the new image
plot(b56r, main="Raster for Lower Teakettle \nBand 56")
Use NEON functions for a way easier way to do what we just did!
# install.packages("devtools")
library(devtools)
# install_github("lwasser/neon-aop-package/neonAOP")
library(neonAOP)
# open band function
b55<-open_band(f,
bandNum=55,
epsg=32611)
plot(b55)
# import several bands
bands<-c(58, 34, 19)
# create a raster stack
RGBStack<-create_stack(f,
bands = bands,
epsg = 32611)
plot(RGBStack) # will plot 3 bands separately
plotRGB(RGBStack,
stretch="lin")
# epsg codes http://spatialreference.org/